/*
    Copyright 2006 Patrik Jonsson, patrik@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise; if not, write to the Free Software
    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
*/

// $Id$

#include "montecarlo.h"
#include <vector>
#include "blitz/array.h"

using namespace std;
using namespace mcrx;

template <typename T>
void test_sampling(T s)
{
  global_random rng;
  seed(1);
  const int n=10;
  array_1 w(n);
  double sum=0;
  for (int i=0;i<n; ++i) {
    w(i)=rng.rnd()*(1-sum);
    sum+=w(i);
  }
  //w= 0.1, 0.25, 0.25, 0.4;
  vector<int> o(n, 0);
  
  s.set_weights(w);
  const int n_sample=1000000;
  for (int i=0; i<n_sample; ++i)
    ++o[s.draw_entry(rng)];

  for (int i=0; i<n; ++i) {
    cout << i << '\t' << w(i) << '\t' << 1.0*o[i]/n_sample << endl;
  }
}

int main(int, char**) 
{
  cout << "Testing cumulative sampling" << endl;
  test_sampling(cumulative_sampling<global_random>());
  cout << "Testing rejection sampling" << endl;
  test_sampling(rejection_sampling<global_random>());
}
